Lack of energy equipartition in homogeneous heated binary granular mixtures 
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We consider the problem of determining the granular temperatures of the components of a homo- 
geneous binary heated mixture of inelastic hard spheres, in the framework of Enskog kinetic theory. 
Equations are derived for the temperatures of each species and their ratio, which is different from 

■ unity, as may be expected since the system is out of equilibrium. We focus on the particular heat- 
ing mechanism where the inelastic energy loss is compensated by an injection through a random 
external force ("stochastic thermostat"). The influence of various parameters and their possible 

I ■ experimental relevance is discussed. 

<n : 

r— i. I. INTRODUCTION 

q ■ ... n 

Experimental and theoretical studies of rapid granular flows |1| have hitherto mostly focused on assemblies of 
identical particles, either freely cooling when the energy loss due to inter-particle collisions is not compensated for, or 
driven in a non-equilibrium stationary state by various energy injection mechanisms. Recently however, interest has 
grown for the more complicated case of polydisperse systems [|, ||, ||, [|, ||, 0, |j| p| [L0|]. Theoretical investigations into 
the homogeneous cooling stage of a binary mixture [|[ ]j| have shown that the two components have different granular 
temperatures (i.e. kinetic energies), even if their cooling rates are equal. Such a result, confirmed by detailed Monte 
Carlo simulations || is also obtained when the system is sheared |], |h]] , heated by the contact with an elastic granular 
gas maintained at fixed temperature T JTl] | , or within the Maxwell model framework p2| . Similarly, a tracer particle 
undergoing inelastic collisions with an equilibrium fluid at temperature T reaches a granular temperature lower than 
/ 0- 

This violation of equipartition in a mixture, although in sharp contrast with the behaviour of molecular gases 
at equilibrium, is not unexpected: the terminology "granular temperature" for the kinetic energy of a granular gas 
has been coined from the equivalence of temperature and kinetic energy in an elastic gas, but does not have any 
thermodynamical status in out-of-equilibrium systems like inelastic granular gases. 

In recent experiments, the granular temperatures have been measured for binary mixtures, both in 3D vibro- 
fluidized granular beds || and in 2D strongly vibrated granular gases . Both studies reported a clear violation of 

■ equipartition with a temperature ratio quite insensitive to the relative densities of the two species. 

£NJ ', The present article aims at providing a simple theoretical framework where the temperature ratio is readily obtained 
£NJ ■ in a non-equilibrium steady state (NESS) . This allows to investigate the influence of many parameters which can be 
difficult to tune experimentally, such as the masses, sizes, densities and inelasticities of the beads. We consider 
analytically a heated binary mixture in the framework of the homogeneous non-linear (Enskog-)Boltzmann equation 
for smooth inelastic hard spheres. Similarly to the case of free cooling described in || and restricting to Gaussian 
velocity distributions, we derive in section II equations for the granular temperatures of the mixture components, 
which are easily solved numerically. The corresponding temperature ratio is in excellent agreement with existing 
numerical work |ll|. In section |lj, we consider the NESS sustained by heating through random kicks ("stochastic 
i-^j ■ thermostat" approach) , a mechanism which has focused some attention recently for one-component (monodisperse) 
systems |lj, [f^, [f7, 18, l|, |2^, ^l], p2|. Although finding an energy injection mechanism of experimental relevance 



is a difficult issue, we expect the approach proposed here to elucidate the basic trends of grain behaviour when varying 
q ' the controlling parameters. Moreover, as will be shown below, the temperature ratio we obtain provides a reasonable 
, zeroth order approximation to compare with the experiments reported in jq, m] . 

X: 



II. KINETIC THEORY 

We consider the model of smooth inelastic hard spheres (IHS) undergoing binary, momentum conserving, inelastic 
collisions, in the framework of the homogeneous non-linear Enskog equation. The system is a mixture of two types of 
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IHS, with masses m\ and m^, diameters a\ and a 2 . Three types of collisions may occur so that the mixture is also 
characterized by three different restitution coefficients: an, a 22 , and a.\ 2 = 021- 

The velocity distributions in the homogeneous state fi(v,t), f 2 (v,t), obey the following kinetic equations: 

d t fi(vi,t) = JMfi, fj] + ?fi (!) 

j 

where the Jij describe the effect of dissipative inter-particle collisions, and T fi represents an external forcing which 
injects energy into the system, allowing it to reach a non-equilibrium steady state. The kernels Jy for collisions 
between a particle of type i and a particle of type j are given, in dimension d, by 



Jij[vi\fi, fj] = Xij<r?j 1 J ' dv 2 J da {a ■ v l2 ) i^-fi{v' x )fj{v 2 ) - fi 



W/iW • (2) 



where the Xij are the pair distribution functions at contact (a priori unknown, but becoming close to 1 in the limit of 
low densities) ; a is a unit vector directed from the center of the particle of type i to the center of particle j [separated 
at contact by CTy = ((7j +aj)/2], and the prime on the integral is a shortcut for J <d(a-vi 2 ). Moreover, V\ 2 = vi — v 2l 
and the pre-collisional velocities v[ and v 2 are given in terms of the post-collisional velocities V\ and v 2 by: 

v[ = vi- + ar. x )(a ■ v 12 )a (3) 
v' 2 =v 2 + + ay 1 )(S : ■ v 12 )a (4) 

where fly — m,i/(mi + to,-), so that momentum is conserved but energy dissipated. 
The partial granular temperatures are defined from the kinetic energies by 

-J-2i(t) = J dv-^-f t (v,t) , (5) 

7ii = J dvfi(v, t) being the number density of particles of type i with a total temperature of the mixture 

T=— ^Vn.T,. (6) 
ni + n 2 

Prom (|l|), the evolution equation for the temperatures reads 

d t Ti = ^ V / dvv 2 J l3 [v\fi, f s ] + J 7 T l , (7) 

J 

where J-Ti describes the effect of the forcing (source) term Tfi. It is possible to integrate over a (see calculations in 
the appendix, and also ||) to obtain, without any further approximation at this stage: 

d-l 



dm = ^i- ^^ 1 ^ J dv 1 dv 2 vf 2 f i (v 1 )f i (v 2 ) 



riid 



d-l 



dv 1 dv 2 f i (v 1 )f j (v 2 ) [fi%{l - a^)vf 2 + 2/^(1 + a t3 )vi 2 {v l2 ■ Vy)} (8) 



with f% = tt^- 1 )/ 2 /T[(d + 3)/2], Vy = fiyVi + fi 3i v 2 , and T the Euler function. 

Since the system reaches a stationary state where the forcing term balances the dissipation due to collisions (the 
forcing and dissipative terms in (j^) generically involve different powers of the temperatures so that its right-hand-side 
admits a "physical" root), we can write dtTi = 0. It is moreover convenient to scale the velocities with the thermal 
velocities Uq,, = y 2Ti / nii , and introduce the functions <&i such that 



fi(v) = — . (9) 

v 0,i \ v 0,iJ 

Equation (|8|) may then be cast into an equation for the rescaled velocity distributions Qf, no further step can however 
be taken without some approximations on the unknown distributions $j . It is convenient to study the deviations of $j 
from the Gaussian $°(c) through an expansion in Sonine polynomials [p3|. In single component heated systems, the 
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deviation from a Gaussian remains small, especially for experimentally relevant values of the restitution coefficient 
jf6], [f8|, |2(], |2l]| . We will here limit our treatment to the Gaussian approximation <E>i(c) = 3>°(c) = ■K d l 2 exp(— c 2 ) 
(lowest order Sonine expansion), ft would of course be possible to go further in a systematic and controlled way, as 
in ||, but we will see by comparison of our approximate analytical calculations with Monte Carlo simulations that, 
at least in the cases we consider, the Gaussian approximation provides reliable results. 

Assuming Gaussian velocity distributions, it is now possible to carry out the remaining integrations in (^); the 



calculations are straightforward and some technical details may be found in the appendix of 
resulting equations for the granular temperatures Ti in the NESS: 

4), 



We only give the 



dT{d/2) 

n , 7r (d-l)/2 



X, 



d-i 2(1 



3/2 
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+ fSL) 
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Mji(l 



2T,; 

rrii 



2Ti 



1/2 



(10) 



These equations still depend on the particular heating mechanism through the term JFT \; once the latter has been 
specified, two equations are obtained for T± and T 2 ] they are easy to implement and solve numerically varying the 
various controlling parameters. 

Before turning to the heating provided by a stochastic thermostat [which amounts to [TTij/rrii — constant], we 
consider three particular limiting cases. 

a. In the tracer limit i.e. m — > 0, without any forcing term, T 2 is imposed and only the equation for T\ is 
considered. As already noted in ||, the result for 7 = T\jT 2 obtained in jl3| 

l + ai2 , , 

7 = — m, T, — (H) 



^(1 



"12 J 



is easily recovered, irrespective of dimension. 

b. Another possibility to obtain a NESS for IHS has been proposed in [ pl| : the temperature T 2 — T of one 
population is imposed, with a corresponding Gaussian velocity distribution, and elastic collisions between particles 
of type 2 as well as for 1 — 2 collisions: 0122 = ol\i = 1. Energy is consequently injected into the inelastic population 
1 (with restitution coefficient an = a < 1), without the need for any other forcing term. In |]TT|| , high precision 
numerical results were obtained for the distribution function $1, and temperature T± (the system is three-dimensional, 
and Xij = l)i by an iterative numerical resolution of the Boltzmann equation. Imposing TT\ = in (|l0|), it is 
straightforward to obtain a third order polynomial equation for 7 = T\jT 2 (this quantity is necessarily smaller than 
1) 
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mi m 2 
— + — + 2 

m 2 mi 



Y = (i-iY 7 + — 



m 1 



m 2 



where 



4ni 



n 2 (l + 0-2/0-1) 



(12) 
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FIG. 1: Comparison of the simulation results found in [fll| (lines) with the solution of equation ( |l2|) (circles), for various values 
of the mass ratio and of the parameter e recalled in (jl^JT 

In Fig. |l|, the solution of equation (|l^) is compared to the results reported in pj| . The agreement is excellent, which 
may be traced back to the analysis of O], showing that the distribution function <!>i is very close to a Gaussian, 
although mathematically different (on the other hand and by definition of the model, $2 is strictly speaking a 
Gaussian). The slight discrepancy obtained at low a for mi = TO2 corresponds to values of the parameters for which 
the deviation of $1 from a Gaussian is stronger. 
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c. Finally, a forcing term Tfi(v) = C£- ■ [vf(v)], which provides an Enskog-Boltzmann equation formally 
equivalent to the free cooling case (see e.g. pjf), leads back to the results of H obtained in this situation: the term 
TT. t is indeed proportional to Tj, so that writing d t Ti = in (|J) yields the same equation for 7 as equating the two 
cooling rates dtTi/Ti when TT i = 0. 



III. THE STOCHASTIC THERMOSTAT 



In this section, we consider the situation of energy supply through random kicks [D| [T(| |T^, |l8|, |l9[ |o[ p2| : 
the particles are submitted between collisions to an uncorrelated white noise (e.g. Gaussian). The equation of motion 
for a particle is then 



dv 

mi — = Fi 
at 



(13) 



where Fi is the force due to inelastic collisions, and (£,ia(t)^jf3(t')) = £,QSijS a pd(t — i'), where Greek indices refer to 
Cartesian coordinates. The associated forcing term in the Enskog equation is 



f(v,t) 



so that TTi = rrii^Q. We do not claim that the forcing term considered here is the most suited to describe vibro- 
fluidized beds, but it mimics an important effect of energy injection by a moving piston: in experiments, particles 
undergoing collisions with the piston (of large mass) gain a velocity that is decorrelated from their masses, so that 
more kinetic energy is injected into the population of large mass. 
The corresponding equation for 7 = T1/T2 reads: 



n 2 \mxj 



(1 



2 _ n i 2 
A 4 21 M12 
n 2 



m 2 
mi 



3/2 



/ n l \ m 2 

2/x 2 i(l + a 12 ) H21 H M12 H 7 

n 2 / V m i 



1/2 



(7-1) 



= X22CT22 



(14) 



The temperature ratio 7 therefore depends in a non-trivial way on the ratios of masses, densities and diameters, 
and also on the inelasticities a%j and pair correlation functions Xij ■ It may be checked that in the limit of vanishing 
inelasticities, 7 — > 1 as it should. Moreover, for mechanically equivalent particles (i.e. mj = m 2 , <J\ = <j 2 and 
an = a 22 — 0112), we should also recover equipartition (7 = 1) irrespective of densities. This is the case in the 
Boltzmann limit (low densities where all pair correlation functions xa ~ * 1)- At arbitrary packing fraction, the 
various approximation for the \ij that may be found in the literature 1 24 , are such that Xij n0 longer depends on 
i and j when <j\ = <J 2 , so that equipartition holds for mechanically equivalent particles. 

Since equation (14) relies on a Gaussian approximation for $i, we have compared our approach to the results of 
Monte Carlo simulations (the so-called DSMC technique [^6)) where the non-linear Boltzmann equation is solved 
numerically for both species. As we solve numerically the homogeneous Boltzmann equation, the phenomena of 
segregation or clustering are explicitly discarded. 

In the following sections, we will study more precisely some cases that could have experimental relevance, and for 
the sake of simplicity, we considered Xij = 1- All the results are given for the three dimensional case; note however 
that for cti = 02 the temperature ratio becomes e?-independent. 



A. Equal inelasticities: ay = a 

We first consider the case of equal restitution coefficients [puj = a), for materials having similar elastic properties. 
The dependence of 7 on the mass and number density ratios, for equal sizes, is shown in Figure ^. Excellent agreement 
is found between DSMC simulations (symbols) and the solution of equation (|l4|). It turns out indeed that the velocity 
distributions measured in Monte Carlo simulations are very close to Gaussians. As may be expected, 7 is a decreasing 
function of m 2 /mi . The density dependence is relatively weak (the temperature ratio slightly increases when n 2 /rii 
increases by an order of magnitude). 

We have also considered two types of beads of the same material, i.e. the same restitution coefficient an = a 22 = ct\ 2 
and same mass density p: the ratio of masses m 2 jm\ is then equal to (o~ 2 /o~i) 3 for three-dimensional beads. Figure 
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FIG. 2: Temperature ratio 7 = T1/T2 as a function of inelasticity, for an = 022 = Q12 = a, and grains of equal radii (01 =02). 
The curves show the solutions of equation ([l4j) whereas the symbols display the results of DSMC simulations. The top three 
curves correspond to a mass ratio 7712 /mi — 2, the three intermediate ones to 7712/7711 = 3 and 7712/7711 = 5 for the three bottom 
curves. For each mass ratio, several density ratios have been considered: 712/711 = 3 (squares and dashed lines), 712/711 = 1 
(pluses and continuous lines), and 712/111 = 1/3 (circles and dot-dashed lines). 

|3| shows a strong influence of the size ratio, for two experimentally relevant values of a: 7 decreases very strongly as 
soon as 02 is two or three times a\. Once again, the number density ratio has a relatively small incidence on 7. It 
is interesting to disentangle the effects of 02/01 and ma/ 'mi, by varying one parameter alone, the other being kept 
constant. It appears that the leading effect in the decrease of 7 observed in Figure || is ascribable to a change in 
mass ratio, and not in size: the results obtained at 01 = 02 varying m%lm\ are close to those reported in Fig. |^, but 
surprisingly give a lower 7 (e.g. the results displayed in Figure || for a = 0.9 and equal densities are 7 = 0.66 and 
0.38 for 02/01 = 2 and 3 respectively, whereas with 01 =02, we obtain 7 = 0.59 for m2/rn\ — 2 3 and 7 = 0.29 for 
m 2 /m 1 — 3 3 ). 




FIG. 3: Temperature ratio from eq (Q) for grains made of the same material, i.e. for 7712/7711 = (0-2/0-1 ) 3 , and otij = a. Filled 
symbols correspond to a = 0.9, open ones to a = 0.7. The squares are for 712/711 = 3, the circles for 712/711 = 1 and the 
diamonds for 7I2/711 = 1/3. 



B. Comparison with experiments 

For glass spheres with size ratio 02/01 = 1.25, Wildman and Parker have measured a temperature ratio 7 = T\jTi 
in the range 0.75-0.8 ||, with a weak dependence on densities (except may be in the limit of large grains predominance 
where 712 ^> rii). Estimating the relevant restitution coefficient to be a ~ 0.9 ||, we obtain from Eq. ( [l4|) 7 ~ 0.9 
(see also Fig. 0), with also a weak dependence on 77,2/711. It is however noteworthy that this weak dependence is 
opposite to that observed experimentally: when the proportion of large grains is increased, we obtain an increase of 
7. For comparison, the temperature ratio obtained for the same parameters in the homogeneous cooling stage |3| is 
7 ~ 0.96 and the authors of || proposed a simplified theory for which 7 is in the range 0.4-0.7. 
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The results obtained by Feitosa and Menon ffi] confirm the very weak influence of n-i / m when the grains (of equal 
size) are made of two different materials. For a mixture glass/ aluminum with mass ratio tt^/toi = 1.09, 7 is measured 
very close to 1, whereas for a more asymmetric mixture of glass and brass with m^/wii ~ 3.6, 7 ~ 0.7. Making use of 
equation (|l4| ) with schematic inelasticities an = ot\2 = 0,21 = 0.85, we obtain 7 = 0.98 for 7712/7711 = 1.09 and 7 = 0.7 
for 7712/7711 = 3.6. In the free cooling regime, the corresponding ratios are 0.99 and 0.82. These results are displayed 
in Fig. I 
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FIG. 4: T1/T2 as a function of mass ratio, for (Xij = 0.85 and o\ = 02, for the stochastic thermostat and for the free cooling 



regime, together with the experimental data of Q. Note that the values at 
as an exact description of the experimental situation. 



0.85 are only schematic and cannot be intended 



Other results are given in Figures |B| and |6| where the partial inelasticities are not taken equal, but are given values 
that we expect to be of experimental relevance: the experimental data of (7) are therefore also reported in Fig. ||(b). 
In Fig. [| the sizes of the particles are taken equal and the mass ratio changes, while Fig. |^ displays the influence of 
the size ratio when the density ratios are fixed. 

The situation reported in J?J corresponds to that of the Figures |^b and |(]b, where the heavier grains are also the 
more dissipative. As may be observed in Fig. |5|b, a variation of ni/772 from 1/3 to 3 leaves 7 roughly unaffected for 
7712/7711 < 3. 

It may be noted that 7 is not bounded from above by 1, and values slightly above 1 are obtained even when 
7712 > tii by conveniently choosing the inelasticities (or, at fixed inelasticities and densities, by conveniently choosing 
the sizes). 7 in nevertheless generically smaller than 1 for 7712 > ttii: the heavier particles have a larger kinetic energy. 
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FIG. 5: (a): T1/T2 as a function of mass ratio, for an = 0.7, ai2 = 0.8, 0^22 = 0.9 and <ti = 02- 

(b): same with "reversed" inelasticities (an = 0.9, 012 = 0.8 and 022 = 0.7), together with the experimental values of |7J. 



IV. CONCLUSION 



We have considered heated binary granular mixtures from the point of view of kinetic theory. As in the free 
cooling case, and in agreement with recent experimental data, the granular temperatures of the components of the 
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FIG. 6: (a): 7 versus size ratio for an equimolar mixture (722 = n{) with an = 0.7, Q12 = 0.8, 022 = 0.9 and different mass 
density ratios pi /pi . 

(b): same with an = 0.9, ai2 = 0.8, 022 = 0.7 and again 712 = n\. 

mixture differ. This finding is not surprising in a non-equilibrium system, where the "temperature" does not have any 
thermodynamical relevance. 

Using a mean-field approach with the assumption of isotropic Gaussian velocity distributions, we have derived an 
equation for the temperature ratio 7 that may be adapted for various kinds of heating mechanisms, and easily solved 
once the controlling parameters have been chosen. In particular, the values obtained within the stochastic thermostat 
framework are compatible with those measured in the experiments reported in || [?j . Even if a quantitative comparison 
with experiments is somehow pointless given the simplicity of our approach, similar trends are observed. For example, 
the heavier particles carry generically more kinetic energy than the lighter ones, the ratio being insensitive to the 
relative number fraction of both species. It also appears that the breakdown of energy equipartition is all the more 
pronounced as the mass ratio is increased, the size ratio playing only a minor role. 

Acknowledgments: We would like to thank T. Biben, R.D. Wildman and CM. Hrenya for communicating their 
results prior to publication. 

APPENDIX 

In this appendix we show how to perform the integrals over cr, in order to obtain equation (|^). We start from the 
identity 



f dvv 2 Jij[v\fi,fj] = Xijvfj 1 J dv 1 dv 2 J da{a ■ Vi 2 )f l {v 1 )f : j{v2){v'( 2 ~ v\) , 



(15) 



with v" = Vi — fJ.ji{l + aij){a ■ V\ 2 )a, i.e. where 



v f -v\ = ^(1 + a l3 ) 2 {S ■ v 12 f - 2/^(1 + aij){a ■ v 12 )($ ■ «i) • 

Using the unit vector C\ 2 = Vi^jvyi., and the known integrals /3„ = J' da(a ■ Cyif 1 (see e.g. flJI), the first term is 
readily computed and yields: 



&Xii<4 Vii(l + Oij) J dv 1 dv 2 f l {v 1 )f J {v 2 )v'l 2 . (16) 
To compute the term containing (<x • vi), we choose one of the unit vectors to be along Vi 2 , and decompose: 



Vl ' V\ 2 ^ I ^ cr ■ V\ 2 ^ ^ I 

vi2 = vi 2 ei, vi — ex+Vi, cr = — e x + a . (17) 



V\1 



V\1 



(cr • Vi) is then written as 



(i>i • v 12 )(S ■ V12) 



cr •«! , 
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and the term tr ■ gives a vanishing contribution in the integral over a for symmetry reasons. We are therefore 
left with 

/ acr 2 = v 12 P3{Vi ■ V12) ■ 

J V 12 

Rearranging terms and writing iq = [ijiViz + fajVi + (ijiV2 — [ijiVvi + Vij, one finally obtains equation (||). 
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